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I. INTRODUCTION 


A. A GENERAL-PURPOSE FINITE ELEMENT PROGRAM 
A general-purpose finite element program should be able to solve a variety of 
problems from the number of disciplines: linear and non-linear, static and dynamic 
problem of elasticity, fluid mechanics, heat transfer, etc. and can solve problems of 
large size involving a variety of elements. 
A general program is going to be voluminous and complex. It is, however, 
desirable that: 
e its logic be easily understood; 
e one or many of its parts be easily modifiable; 


e it os possibilities to tailor its facilities for the solution particular classes of 
problems. 


The program should have a modular structure, with the modules made as 
independent from one another as 1s practicable. The following modular operations are 
recognized: 

1. Problem Definition (data base): 

- node coordinates and element connectivities; 
- nodal element properties; 


- boundary conditions. 


tJ 


. Element Computations: 
- integration points and associated weights; 
- interpolation functions and their derivatives; 
- Jacobian matrix, inverse, and determinant; 


- element matrices and vectors: [k], [ml], {f}, etc. 


re] 


. Assembly Operations: 


- assemblage of master matrices and vectors, [K], [Ml], etc. 


es 


. Solution: 
- factorization of master matrices and solution of equations. 
5. Result: 
- output of nodal variables and other calculated quantities: gradients, 


Teactions cic. 


Subroutines implementing the various operations described above are 
contained in all finite element codes. The flow of information between these 
operations is problem dependent; linear, non-linear, static, and dynamic problems all 


require logic of their own. 


B. MEF PROGRAM 

A program of medium complexity, called MEF, implementing the techniques of 
general-purpose program that can solve a large variety of boundary value problems of 
mathematical physics. It is written in FORTRAN IV and can be easily adaptable to 
various computers. 

The main program controls the flow of all information through the functional 
blocks by transferring control to a subroutine called BLNNNWN when the block calling 
card NNNN is encountered in the input file. The subroutine BLNNNN then performs 
preliminary functions such as logical unit identification, and reading of control 
parameters for the creation of various files and tables. The subroutine then calls 
subroutine EXNNNN. In all cases, subroutine BLNNNWN provides appropriate default 
parameters which will be overridden by user values if specified. Subroutine EXNNNN 
then performs the major operations of the block by calling on the needed subroutines 
in the MEF library. The above protocol holds for all blocks except STOP, COMT, 
and IMAG. All the functions of COMT and IMAG are performed by subroutine 
BLNNNN, and the function of block STOP is performed by the main program. 


The executable functional blocks contained in the MEF are: 


BEOCK FUNCTION 

SOLR ~~ Assemblage of distributed load. 

LINM — Solution of linear problem with global matrix in core. 
LIND — Solution of linear problem with global matrix out of core. 
NLIN — Solution of stationary non-linear problem. 

TEMP — Solution of unsteady problem (linear or non-linear). 


VALP Eigenvalues and eigenvectors. 


The various blocks designed for execution of the various computations have 
similar structures since they have to: 
® construct element and load matrices: 
e assemble global matrices and vectors; 


e factorize and solve the system of equations; 
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® output the results. 

Using the constructed elements and load matrices, the subroutine element library 
ELEMLB 1s called. This library contains subroutines that define the individual element 
types. The ELEMO3 subroutine, a thirty-two node, three dimensional isoparametic 
element was developed and added to the element hbrary. This new element allowed the 
solution of linear elastic structures composed of homogeneous and isotropic materials. 

Multiple sample problems were developed to fully exercise use of this new 
element. An indepth investigation was then conducted to determine the limit of 
computational ability using the newly defined element to represent physical 


phenomenons. 
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II. DESCRIPTION OF THE COMPUTER PROGRAM 


A. REFERENCE ELEMENT 

To simplify the analytical expression for elements of complex shapes an element 
of reference is introduced. Such an element is defined in an abstract non-dimensional 
space with a very simple geometrical shape. The geometry of the reference element is 
then mapped into the geometry of the real element using geometrical transformation 
expression. 

A thirty-two node, three dimensional cubic element was introduced as the 
reference element. The element has eight coner nodes and twenty four mid-side nodes 
dividing each edge in three equal parts as shown in Figure 2.1. Using this reference 
element we created the fundamental matrices and vectors in subroutine called 
ELEM03, NI03, D0O3 and BO3 to be used in element library subroutine ELEMLB of 
tac MEF program. 





Figure 2:1 Reference Element. 


B. SUBROUTINE ELEMO03 

Multi-purpose program MEF 1s organized to contain a library of one, two, and 
three-dimensional elements for the solution of problems from a wide variety of 
disciplines. Problems from the mechanics of solids and fluids, heat transfer, etc. have 
been solved. For each element type mn, subroutine ELEMnn controls the computation 
of all matrices and vectors. 

In element type 03, subroutine ELEM03 is organized to create the fundamental 
matrices and vectors that must be numerically integrated using methods of Numerical 
Integration described in and the computations can be carried out by the following 
steps. 

1. Operation common to all element of the same type: 

- compute the weight w, and the coordinate of integration points; 
- construct the functions N (interpolation functions), the function N 
(geometrical interpolation functions) and their derivatives 
with respect to 6, n, § at the points of integration. 
2. Operation for the computation of matrix [k] of each element: 
- initialize the matrix [k]; 
- for each point of integration 6, : 
¢ construct the Jacobian matrix [J] from the derivatives with respect 
to &, n, § of function N and the nodal coordinates of the element; 
® construct the inverse of [J] and its determinant; 
® construct the derivatives of functions N with respect to x, y, z 
starting from the derivatives with respect of &, n, 6; 
® construct the matrices [B] and [D]; 
¢ accumulate into [k] the values of [B]'[D][B]det(J)w, calculated for 
each integration point. 
3. Operation required to compute mass matrix [m]: 
- initialize [m]} 
- for each integration point §, : 
¢ compute Jacobian matrix and its determinant; 
¢ accumulate the values of {N} < N> det(J)w, into matrix [m]. 
4. Operation required to compute consistent load vectors {f}: 
- initialize {f}; 
- for each integration point 6, : 
¢ compute the Jacobian matrix and its determinant; 
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¢ accumulate the values of {N}f,, det(J)w, into {f}. 
5. Operation required to compute the residue vecter {r}: 
- using the value of {f} from 4; 
- for each integration point ¢, : 
® compute matrices [B], [D], [J] as in 2 above; 
¢ accumulate the product: {f} - [B]'[D][B]{u,} w,det(J) into {r}. 
6. Operation required to compute gradients du at points of integral: 
- for each integration point &, : 
® construct matrix [B] as in 2 above; 
¢ compute and print gradient {du} = [BJ]{u,}. 
The subroutine ELEM0O3 executes one operation at a time depending on the 
value of ICODE. Control variable ICODE specifies which element operation is 


desired, and expression of this variable as follows: 


ICODE = | initialization of the characteristic parameters of an element 
(number of nodes, number of degrees of liberty). 

ICODE = 2 operations required by a given referance element which are 
independent of the real geometry; construction of 
interpolation function N and their derivative with 
respect to & at the points of integration. 

ICODE = 3 construction of matrix [k] in array VKE. 

ICODE = 4 construction of tangent matmx needed for non-linear 
problems in array VKE. 

ICODE = 5 construction of mass matrix [m] in array VKE. 

ICODE = 6 computation of residual vector {r} in array VFE. 

ICODE = 7 computation load vector {f} in array VFE. 

ICODE = 8 computation and printing of gradients {du}. 

C. CODING. 


1. Evaluate coordinates, weights, functions N and their derivatives 


Integration formular for numerical integration is the following form. 


Lng 
es 2 ep (eqn 2.1) 
where a 
- Gj are the coordinates of integration point ‘I’ in §, n, § system 


coordinate corresponding to weight whys 
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- W, are the weights corresponding to integration point number, 
- IPG are the total number of integration points. 

A choice of 2, 3 or 4 integration points by dimension can be made in 
subroutine GAUSS [Ref. I:p. 265], giving respectively 8, 27 and 64 integration points, 
the weights corresponding to the integration points and their coordinates. They are in 
the array called IPG, VCPG and VKPG respectively. 

Subroutine NIQ3 create the array VNI that contains the shape function N; 


and their derivatives with respect to §; ( §, n, § system coordinate ) as Figure 2.2. 


For each reference element 


For each integration 
“point ¢, 





Figure 2.2 Block Diagram for the Shape Function and their Derivative. 


2. Computation of stiffness matrix [k] of each element. 
The explicit equation of element stiffness matrix 1s following: 
le ee: 
aoe [B]'[D][B]det(J)dg dn d& (eqn 2.2) 
where 
- [B] is the linear strain matrix; 
- [By]! is the transpose matrix of [B]; 
- [D] is the matrix of elastic constants for an isotropic material. 
The stiffness matrix that has the block diagram as Figure 2.3 is the integration 
of product of the linear strain transpose matrix, the element property matrix and the 


linear strain matrix over the volume of the reference element. 
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For each integration 
point € 










Subroutine JACOB 
(JJ, (3), det(J) 






Subroutine DNIDX 
ON; ON; ON; 
Ox Gy az 


aeme:). 


subroutine D03 
[D] 












Subroutine BTDB 
[k] = [B]{D][B]det(J) 





Figure 2.3 Block Digram of Computation of Stiffness Matrix. 


In coding, we employ numerical integration formular having the following 


generic form: 


where 


ipg 
| = 2 (S))] (eqn 2.3) 


- IPG are the total number of integration points; 

- W; are the weighting coefficients corresponding to each integration 
point; 

- §; are the coordinate of the IPG integration points; 


- {k°] is the stiffness matrix at each integration point as shown in 


equation 2.4. 


IS 


(k°] = [B]{D][B]det(J) (eqn 2.4) 


Subroutine JACOB [Ref. I:p. 63], compute the Jacobian matrix, its inverse 
and its determinant. The Jacobian matrix as Figure 2.4 is obtained as the product of 
two matrices, one containing the derivatives of the geometrical transformation 
functions with respect to the space of the reference element, and the other containing 


the real coordinates of the geometrical nodes of the element. 


<xy2> = <NG)= ey (a (eqn 2.5) 


{Xn}, {Yn}, {Z,} being the geometrical node coordinates. The Jacobian matrix is figure 
24 


Ox Gy 2 
a GE 
i= 9S Se 
on on an 
Ox dy @2 
a ae a 


Figure 2.4 The Jacobian Matrix. 


The inverse of Jacobian matrix ay! as Figure 2.5 and its determinant are 
computed by the method discribed on [Ref. l:p. 44], 

Subroutine DNIDX [Ref: l:p. 64], computes the derivatives of the shape 
function NV, with respect to the coordinate system of the real element using the product 
One Figure. 2-0. 

Subroutine BO3 creates the matrix as Figure 2.7 that contains the strain 
components at all direction of each node in the element using output components of 
subroutine DNIDX and rearrange the new array called VBE. 

Subroutine D0O3 compute the stress-strain matrix (VDE) as Figure 2.8 


{Ref. 2:p. 110] for isotropic materials ( E = Young’s Modulus, v = Pisson’s Ratio ). 
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Subroutine BTDB construct the stiffness matmx by adding the product of the 
transpose matrix VBE, the stress-strain matmx VDE and the matrix VBE of every 


integration points. 
[k] = [B]'{D][B]w,det(J) (eqn 2.6) 


3. Computation The Mass Matrix [m]. 


The element mass matrix has the block diagram as Figure 2.9 and the explicit 


equation as following: 


els) 
[m] = J j j (N]{N]det(J) do dn df (eqn 2.7) 


and numerical integration form Is: 





Figure 2.8 The Array VDE: 


ipg 
[m] = + w;[N]'[N]det(J) (eqn 2.8) 
i=l 


- [N] is the shape function matrix Figure 2.10 that was created by 
subroutine NIO3; 
- [N]' is the transpose matrix of [N]: 
- W; are the weighting coefficients corresponding to each integration 
point; 
Obviously look at the product of matrices [N]‘and[N] as Figure 2.11 is a 
symmetric matrix. By convention, the mass matrix [m] contains upper half components 


and diagonal components of this matrix. 


18 


For each integration 
point Sr 

Subroutine JACOB 

: ale (JJ, det(J) 


Accumulate the values 







[m] = w;[N]'[N]det(J) 


Figure 2.9 Block Diagram for the Mass Matrix [ml]. 





Figure 2.10 The Shape Function Matrix [N]. 


4. Computation of consistent load vector {f}. 


The explicit formular for the load vector {f} is: 


piu exe 
= J J J [N] {f,}det(J) dO dn dé (eqn 2.9) 
VZ 


has the block diagram as Figure 2.12 and numerical integration form is: 


ipg tx 
{f} = ) wi{N]{f vy del) (eqn 2.10) 
i=] free , 


IU, 
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Figure 2.11 The Product of Matrices [N]° and [N]. 


eee fy Ge are the force per unit volume in direction x, y, z. 
5. Computation the residue array. 
The element residual array has the block diagram as Figure 2.13 and defined 


by: 
(yeahs Key, (eqn 2.11) 


where 
- {r} is the residue vector; 
- {f} 1s the element load vector; 
- [k] is the element stiffness matrix; 
- {u,} is the nodal values vector. 
6. Computation of the gradients and stresses at the points of integration. 
For each integration point: 
- compute strain as Figure 2.14. [Ref. 3:p. 31] 


or compacted form 


{c} = [Bl{u,) (eqn 2.12) 
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For each integration 
‘point &. 


Subroutine JACOB 
fad, oly det(J) 
















Accumulate the values 
{f} = {N}f, detJ)w, 


Figure 2.12 Block Diagram for Consistent Load Vector [f]. 


For each integration 





point &§- 





Subroutine JACOB 
(JJ, [J], det(J) 





Subroutine D03 
Subroutine 503 


(D], [3B] 





Accumulate the values 
fr} = {f} - [BI'{DI[B]{uy} w,det(J) 


Figure 2.13 Block Diagram for the Residue Array. 


- compute stress as Figure 2.15. [Ref. 3:p. 30] 


or compacted form 
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Figure 2.14 The Strain-Node Value Relation. 





Figure 2.15 The Stress-Strain Relation. 


{o} = [D]t&} (eqn 2.13) 


- compute coordinate of each integration points on the real element as 
figure 2.16. 


or compacted form 
GS ier (eqn 2.14) 


Using X,, Yq, Z, »where N = N(§;, nj, 6) and Si. Ny G; are the gauss points, 
the coordinate of thirty two nodes on the real element that existed in the array 
VCORE. All the integration points are evaluated in the reference element. Using the 
product of the transfer matrix (the shape function matrix N) and the array VCORE, we 


can evaluate the coordinate of the integration points on the real element. 


tN 
tr 





Figure 2.16 Evaluation of the Coordinate of the Integration Points on the Real Element. 


Note that subroutine ELEM0O3 executes one operation at a time depending on 
the value of ICODE. For preserving the memory locations, some of the operations 
have been performed using the same array name as VKE in both the stiffness matnx 
[k] and the mass matrix [m]. The same thing has been done with array VFE used for 


the residual vector {r} and the load vector {f}. 


Ill. THE SOLUTION OF A Sct too AND DISCUSSION OF 


In this chapter the preparation of input data for the computer program are 
described. A simple problem was developed and the investigation was conducted to 


determine the ability of this cubic element. 


A. ENTRY AND EXECUTION FUNCTIONAL BLOCKS 

MEF has specialized functional blocks for the entry, verification and 
Organization of the data required to define a problem. Block COOR reads the nodal 
coordinates and the number of degrees of freedom of each node, it also provides 
automatic node generation. Block COND reads the boundary condition. Block PREL 
reads element properties if required for element type being used. Block SOLC reads 
the concentrated loads. Block ELEM reads the element connectivities; it also reads 
element group information when more than one element type is used, if elements have 
different properties. This block provides automatic element generation. 

Other function blocks of MEF for the execution of particular finite element 
computations use the data base constructed by entry blocks and augment it by their 
results. Block LINM assembles and solves a linear system of equations residing 
in-core. Block LIND is similar to the block LINM but the system of equations resides 
out-of-core in a mass storage device. And block STOP terminates execution of the 
problem. 

MEF provides various levels of output. The quantity of output desired from a 
given block is controlled by a parameter on the block calling card, described in detial 
[Ref: l:pp. 440-447], which ranges from 0 (the assumed value) to 4. The default value 
provides all the information needed to verify the input stream and obtain the desired 


answers while the value | thru 4 provid? various level of verbosity. 


B. STRUCTURE MODELING 

The first step in applying a finite element solution to the problem 1s the selection 
of an appropriate mesh. In many case an extrapolation of the results will be required 
and hence more than one mesh will have to be selected. The mesh size solution does 
not follow any predetermined rules and will, in general, depend on the nature of the 


problem and the judgement of the analyst. An abitrary numbering convention 1s 
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adopted for the nodal points of the entire structure (in distinction with the numbering 
convention for the element nodes shown in Figure 3.1. A second numbering scheme is 
also needed for the elements of a mesh. 

In order to show the function of the program, as well as some observations 
affecting the use of it, the solution of a simple problem is presented in this section. 
The problem selected is that usually presented in classical texts of strength of materials 
as a cantilever beam of uniform cross section subjected to a concentrated load at the 
end of the beam. The model used for this problem is shown in Figure 3.1. Nodes I, 2, 
3, 4, 5, 6, 7, 8, 9, 10, 11 and 12 (here arbitrarily numbered) are constrained to zero 
displacement in all directions. The concentrated load at the end of the cantilever beam 
is considered as the consistent load and are shown in Piptire o.2: 

Four different meshes which each mesh has the thickness 9, 0.9, 0.09, 0.009 
inches and has the concentrated loads 1200, 1.2, 1.2e-3 and 1.2e-6 pounds were 
employed in order to show the variation of results with mesh size. Numerical results 


for the maximum displacement at the tip of the beam are shown in Tables I and II. 





Figure 3.1 The Model used for the Problem. 
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Figure 3.2 The consistent load at the end of the beam. 


C. DISCUSSION OF RESULTS 

Elementary beam theory gives a displacement of 0.0631 inches, however this 
beam theory neglects shear deflection and three dimensional elasticity does not do so. 
Observation of the results presented here reveal that mesh refinement lead to improved 
results which eventually will converge to a certain value. The thirty two nodal points 
brick is subject to numerical bad conditionary when used with an adverse slenderness 
ratio. This ratio being defined as the ratio of the maximum element dimension over 
the minimum dimension. For example tn the numerical examples treated the 
slenderness ratio for the one and eight elements representation of the four beams 
analysis are shown in the Tables I and II. Results are acceptable for the first three 
beams in each case. However a computation of approximate values of zero, first 
invariant of the element stiffness as well as condition number are produced to 
investigate the results. 

As the slenderness ratio increases the numerical of the conditioning of the 
stiffness matrices become so bad that no significant digits can be expected out of the 
solution for displacements. The same conclusion is arrived with a mesh of eight 
elements. In such a case the slenderness ratio varies from 1.6 to 1666.6. For the 
thickness of 0.09 the slenderness ratio was adequate. To reduce the ratio, more 
element must be used and tests with 20 and 40 elements were conducted. The twenty 
elements mesh gives of 0.6 to 666.6. In the use of a mesh with fourty elements the 


largest dimension become 6.0 and must therefore be reduced to 3.0 as 3.0 is the value 
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of 120/40. Runs with such a mesh were then performed to confirm our results. On the 
basis of all these runs and the value of the condition number of these matrices it seems 
that a slenderness ratio of approximately 150 could still be used to give results with 6 
significant digits in the answers, however further examples must be treated before a 
conclusion could be reached. 

These are observations essentially made on the results of a simple problem. Thus 
no firm rules regarding the use of the newly element can be established and further 


experimention with different problems has to be conducted for this purpose. 
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TABLE I 
THE DISPLACEMENT FOR THE ONE ELEMENT 


Slenderness 


Thickness Ritia. - 


@0.54e-5 7.27el11 3.72e10 1.18e06 


-0.0581 


-0.16e-4 7.11e12 1.64e11 2.22209 “0.0515 


0. 74e-4 7.11e13 1.64e12 7.42e<2 2.22e13 


0. 74e-3 0.56e-6 7.11e14 1.64e13 7.42e23 4.59e15 





TABLE I] 
THE DISPLACEMENT FOR THE EIGHT ELEMENTS 


Slenderness Conditton 


Nie Ber | Displacement 


Ratio 


-0.0625 


0. 33e-6 


0. 94e-6 


0.93e-5 


0.93e-4 


0.12e-4 


0.97e-4 


3.21410 


9.12e10 


8.89eil 


8.88e12 
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4.66e09 


2.06e10 


2.06e11 


2.06e12 


1. 86e06 


1.60e04 


1.65e01 


2.-54e03 


1.29e06 


1.29e10 


1.32e14 


-0.0614 


-0.0597 





APPENDIX 
ELEM03 SUBPROGRAM LISTING 


It includes four 


The listing of subprogram ELEMO3 is provided below. 


subroutines: NIO3, DO3, BO3 and BTDB respectively. 
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s----— CHARACTERISTIC DIMENSIONS OF THE ELEMENT 


VI(NDIM**2) , VJ1(NDIM**2) 
INEL* IPG) , IPGKED(NDIM) 


DE1( IMATD**2) 


------- DIMENSION OF MATRIX D, NUMBER OF G. P. 


a 05/0 soUG/ SRADN/ 57 207 795160823027 
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tao co cceeeecoose FUNCTION 70 BE EXECUTED 
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